function fig1_11

nr=50;
rMax=10;
alpha2=0.1; z1=1; z2=-2; alpha1=-alpha2*z2/z1;
k=sqrt(alpha1*z1^2+alpha2*z2^2);
lambda=0.5*(alpha1*z1^3+alpha2*z2^3);
ep=0.1;

r=linspace(1,rMax,nr);
gamma=lambda*exp(2*k)/(2*k)/(1+k)^2;
a=gamma*( (k-1)*exp(2*k)*expint(3*k) + (k+1)*expint(k) )/(k*(1+k));

for ir=1:nr
	ar=k*r(ir);
	phi0(ir)=exp(k-ar)/r(ir)/(1+k);
        phi1(ir) = a*exp(-ar)/r(ir) + gamma*( exp(ar)*expint(3*ar) - exp(-ar)*expint(ar) )/(ar);
end

phi=ep*(phi0+ep*phi1);

% note rM>rMax
rM=20;

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.2,'MR',-0.01,'ML',0.058,'P',0.04)

%  calculate solution using a MATLAB routine
sol0 = bvpinit(linspace(1,rM,1000),@mat4init);
sol = bvp4c(@rhs,@bcs,sol0)

%  plot 
hold on
xint = linspace(1,rM);
Sxint = deval(sol,xint);

box on
plot(r,ep*phi0,'r','Linewidth',1.2)
plot(xint,Sxint(1,:),'--b','Linewidth',1.4)
box on
grid on
axis([1 rMax 0 0.06])
ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$r$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('Potential','Interpreter','latex','FontSize',14)
legend(' $\, \varepsilon\phi_0 \,$','$\, \phi \,$','Location','NorthEast','FontSize',13,'Interpreter','latex')

hold off
%exportgraphics(gcf,'/Users/mark/Desktop/pb.pdf')



%  define RHS of equations
function dy=rhs(x,y)
dy=zeros(2,1);
dy(1) = y(2);
dy(2) = p(y(1)) - 2*y(2)/x;

%  define BCs
function res = bcs(ya,yb)
ep=0.1;
res = [  ya(2) + ep
         yb(1)  ];

function yinit = mat4init(x)
ep=0.1; 
alpha2=0.1; z1=1; z2=-2; alpha1=-alpha2*z2/z1;
k=sqrt(alpha1*z1^2+alpha2*z2^2);
c = ep*exp(k*(1-x))/((1+k)*x);
yinit = c* [  1
          -k-1/x ];


function g=p(y)
alpha2=0.1; z1=1; z2=-2; alpha1=-alpha2*z2/z1;
g= - alpha1*z1*exp(-z1*y) - alpha2*z2*exp(-z2*y) ;
































